home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / exampleCode / games / IndiZone / blix / blixvect.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  13.6 KB  |  616 lines
  1. /*
  2.  * Copyright (C) 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*_________________________________________________________________________
  18.  |
  19.  | blixvect.c - functions to support operations on vectors and matrices.
  20.  |
  21.  | Lots of original code from David Ciemiewicz, Mark Grossman, Henry Moreton,
  22.  | Paul Haeberli, and Gavin Bell.
  23.  |
  24.  | Frans van Hoesel added more specific prototypes,
  25.  | and some functions (not all tested)
  26.  |
  27. */
  28.  
  29. #include <malloc.h>
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <bstring.h>
  33. #include "blixvect.h"
  34. #include "mymath.h"
  35.  
  36.  
  37. Matrix idmatrix =
  38. {
  39.     { 1.0, 0.0, 0.0, 0.0,},
  40.     { 0.0, 1.0, 0.0, 0.0,},
  41.     { 0.0, 0.0, 1.0, 0.0,},
  42.     { 0.0, 0.0, 0.0, 1.0,},
  43. };
  44.  
  45. float x_axis[] = { 1.0, 0.0, 0.0 };
  46. float y_axis[] = { 0.0, 1.0, 0.0 };
  47. float z_axis[] = { 0.0, 0.0, 1.0 };
  48.  
  49.  
  50.  
  51. float *vnew(void) {
  52.  
  53.     return ( (float *) malloc(3 * sizeof(float)));
  54. }
  55.  
  56. float *vclone(const float v[3]) {
  57.     register float *c;
  58.  
  59.     c = vnew();
  60.     vcopy(v, c);
  61.     return c;
  62. }
  63.  
  64. void vcopy(const float v1[3], float v2[3]) {
  65.     v2[X] = v1[X];
  66.     v2[Y] = v1[Y];
  67.     v2[Z] = v1[Z];
  68. }
  69.  
  70. void vcopy4(const float v1[4], float v2[4]) {
  71.     v2[X] = v1[X];
  72.     v2[Y] = v1[Y];
  73.     v2[Z] = v1[Z];
  74.     v2[W] = v1[W];
  75. }
  76.  
  77. void vprint(const float v[3]) {
  78.     printf("x: %f y: %f z: %f\n",v[X],v[Y],v[Z]);
  79. }
  80.  
  81. void vprint4(const float v[4]) {
  82.     printf("x: %f y: %f z: %f w: %f\n",v[X],v[Y],v[Z],v[W]);
  83. }
  84.  
  85. void vset(float v[3], const float x, const float y, const float z) {
  86.     v[X] = x;
  87.     v[Y] = y;
  88.     v[Z] = z;
  89. }
  90.  
  91. void vset4(float v[3], const float x, const float y, const float z,
  92.         const float w) {
  93.     v[X] = x;
  94.     v[Y] = y;
  95.     v[Z] = z;
  96.     v[W] = w;
  97. }
  98.  
  99. void vzero(float v[3]) {
  100.     v[X] = 0.0;
  101.     v[Y] = 0.0;
  102.     v[Z] = 0.0;
  103. }
  104.  
  105. void vzero4(float v[4]) {
  106.     v[X] = 0.0;
  107.     v[Y] = 0.0;
  108.     v[Z] = 0.0;
  109.     v[W] = 0.0;
  110. }
  111.  
  112. void vnormal(float v[3]) {
  113.     float l;
  114.     l =  sqrtf(v[X] * v[X] + v[Y] * v[Y] + v[Z] * v[Z]);
  115.     v[X] /= l;
  116.     v[Y] /= l;
  117.     v[Z] /= l;
  118. }
  119.  
  120. void vnormal2(const float v[3], float dst[3]) {
  121.     float l;
  122.     l = sqrtf(v[X] * v[X] + v[Y] * v[Y] + v[Z] * v[Z]);
  123.     dst[0] = v[0] / l;
  124.     dst[1] = v[1] / l;
  125.     dst[2] = v[2] / l;
  126. }
  127.  
  128. float vlength(const float v[3]) {
  129.     return sqrtf(v[X] * v[X] + v[Y] * v[Y] + v[Z] * v[Z]);
  130. }
  131.  
  132. void vscalar(float v[3], const float mul) {
  133.     v[X] *= mul;
  134.     v[Y] *= mul;
  135.     v[Z] *= mul;
  136. }
  137.  
  138. void vmult(const float src1[3], const float src2[3], float dst[3]) {
  139.  
  140.     dst[X] = src1[X] * src2[X];
  141.     dst[Y] = src1[Y] * src2[Y];
  142.     dst[Z] = src1[Z] * src2[Z];
  143. }
  144.  
  145. void vadd(const float src1[3], const float src2[3], float dst[3]) {
  146.     dst[X] = src1[X] + src2[X];
  147.     dst[Y] = src1[Y] + src2[Y];
  148.     dst[Z] = src1[Z] + src2[Z];
  149. }
  150.  
  151. void vsub(const float src1[3], const float src2[3], float dst[3]) {
  152.     dst[X] = src1[X] - src2[X];
  153.     dst[Y] = src1[Y] - src2[Y];
  154.     dst[Z] = src1[Z] - src2[Z];
  155. }
  156.  
  157. void vhalf(const float v1[3], const float v2[3], float half[3]) {
  158.     float len;
  159.  
  160.     vadd(v2,v1,half);
  161.     len = vlength(half);
  162.     if(len > 0.0001)
  163.     vscalar(half,1.0/len);
  164.     else
  165.     vcopy(v1, half);
  166. }
  167.  
  168. float vdot(const float v1[3], const float v2[3]) {
  169.     return v1[X]*v2[X] + v1[Y]*v2[Y] + v1[Z]*v2[Z];
  170. }
  171.  
  172. void vcross(const float v1[3], const float v2[3], float cross[3]) {
  173.     float t0, t1;
  174.  
  175.     t0 = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
  176.     t1 = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
  177.     cross[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
  178.     cross[X] = t0;
  179.     cross[Y] = t1;
  180. }
  181.  
  182. void vdirection(const float v1[3], const float v2[3], float dir[3]) {
  183.     float t0, t1, t2;
  184.     float l;
  185.     
  186.     t0 = v2[X] - v1[X];
  187.     t1 = v2[Y] - v1[Y];
  188.     t2 = v2[Z] - v1[Z];
  189.     l = t0*t0 + t1*t1 + t2*t2;
  190.     if (l <= 1e-9 ) {
  191.     /* create a direction cross to v1 */
  192.     t0 = 0;
  193.     t1 = -v1[Z];
  194.     t2 = v1[Y];
  195.     l = t1*t1 + t2*t2;
  196.     }
  197.     /* normalize dir */
  198.     l = sqrtf(l);
  199.     dir[X] = t0 / l;
  200.     dir[Y] = t1 / l;
  201.     dir[Z] = t2 / l;
  202. }
  203.  
  204. void vreflect(const float in[3], const float mirror[3], float out[3]) {
  205.     float temp[3];
  206.  
  207.     vcopy(mirror, temp);
  208.     vscalar(temp,vdot(mirror,in));
  209.     vsub(temp,in,out);
  210.     vadd(temp,out,out);
  211. }
  212.  
  213. int veq(const float v1[3], const float v2[3]) {
  214.     if (v1[0] == v2[0] && v1[1] == v2[1] && v1[2] == v2[2])
  215.     return 1;
  216.     else
  217.     return 0;
  218. }
  219.  
  220. /*
  221. * polar coordinates are encoded as eighter theta, phi, or as
  222. * theta, phi, r. The case with only theta and phi is asumed to have
  223. * a radius of 1.0
  224. *
  225. */
  226.  
  227. void frompolar(const float p[2], float v[3]) {
  228.     register float spt;
  229.  
  230.     spt = sinf(p[THETA]);
  231.     v[X] = spt * cosf(p[PHI]);
  232.     v[Y] = spt * sinf(p[PHI]);
  233.     v[Z] = cosf(p[THETA]) ;
  234. }
  235.     
  236. void frompolar3(const float p[3], float v[3]) {
  237.     register float spt;
  238.  
  239.     spt = sinf(p[THETA]);
  240.     v[X] = p[R] * spt * cosf(p[PHI]);
  241.     v[Y] = p[R] * spt * sinf(p[PHI]);
  242.     v[Z] = p[R] * cosf(p[THETA]) ;
  243. }
  244.  
  245. void topolar(const float v[3], float p[2]) {
  246.     register float t;
  247.  
  248.     p[THETA] = acosf(v[Z] );
  249.     t = sqrtf(1 - SQR(v[Z]));
  250.     if (t < 1e-5)
  251.     p[PHI] = 0.0;
  252.     else
  253.     p[PHI] = atan2f(v[Y]/t, v[X]/t);
  254. }
  255.  
  256. void topolar3(const float v[3], float p[3]) {
  257.     register float r;
  258.  
  259.     r = sqrtf(SQR(v[X]) + SQR(v[Y]) + SQR(v[Z]));
  260.     p[R] = r;
  261.     if (r < 1e-5) {
  262.     p[THETA] = 0.0;
  263.     p[PHI] = 0.0;
  264.     } else {
  265.     p[THETA] = acosf(v[Z] / r);
  266.     /* re-use register variable r */
  267.     r = sqrtf(SQR(r) - SQR(v[Z]));
  268.     if (r < 1e-5)
  269.         p[PHI] = 0.0;
  270.     else
  271.         p[PHI] = atan2f(v[Y]/r, v[X]/r);
  272.     }
  273. }
  274.  
  275. float vdistance(const float v1[3], const float v2[3]) {
  276.     
  277.     return sqrtf(SQR(v1[X]-v2[X]) + SQR(v1[Y]-v2[Y]) + SQR(v1[Z]-v2[Z]));
  278. }
  279.  
  280. float vdistance2(const float v1[3], const float v2[3]) {
  281.     
  282.     return SQR(v1[X]-v2[X]) + SQR(v1[Y]-v2[Y]) + SQR(v1[Z]-v2[Z]);
  283. }
  284.  
  285. float vdistance_angle(const float p1[2], const float p2[2]) {
  286.  
  287.     return acos(cos(p1[THETA])*cos(p2[THETA]) +
  288.         sin(p1[THETA])*sin(p2[THETA])*cos(p1[PHI]-p2[PHI]));
  289. }
  290.  
  291. void mprint(const Matrix m) {
  292.     register int row;
  293.  
  294.     printf("Matrix: \n");
  295.     for (row=X; row <=W; row++)
  296.     printf("%12f %12f %12f %12f\n", m[row][X], m[row][Y],
  297.         m[row][Z], m[row][W]);
  298. }
  299.  
  300. void mmultmatrix(const Matrix m1, const Matrix m2, Matrix prod) {
  301.     register int row, col;
  302.     Matrix temp;
  303.  
  304.     for(row=X ; row<=W; row++) 
  305.     for(col=X; col<=W; col++)
  306.         temp[row][col] = m1[row][X] * m2[X][col]
  307.            + m1[row][Y] * m2[Y][col]
  308.            + m1[row][Z] * m2[Z][col]
  309.            + m1[row][W] * m2[W][col];
  310.     for(row=X ; row<=W; row++) 
  311.     for(col=X ; col<=W ; col++)
  312.         prod[row][col] = temp[row][col];
  313. }
  314.  
  315. void vtransform(const float v[3], const Matrix mat, float vt[3]) {
  316.     float tx, ty;
  317.  
  318.     tx = v[X]*mat[X][X] + v[Y]*mat[Y][X] + v[Z]*mat[Z][X] + mat[W][X];
  319.     ty = v[X]*mat[X][Y] + v[Y]*mat[Y][Y] + v[Z]*mat[Z][Y] + mat[W][Y];
  320.     vt[Z] = v[X]*mat[X][Z] + v[Y]*mat[Y][Z] + v[Z]*mat[Z][Z] + mat[W][Z];
  321.     vt[X] = tx;
  322.     vt[Y] = ty;
  323. }
  324.  
  325. void vtransform4(const float v[4], const Matrix mat, float vt[4]) {
  326.     float t[4];
  327.  
  328.     t[X] = v[X]*mat[X][X] + v[Y]*mat[Y][X] + v[Z]*mat[Z][X] + mat[W][X];
  329.     t[Y] = v[X]*mat[X][Y] + v[Y]*mat[Y][Y] + v[Z]*mat[Z][Y] + mat[W][Y];
  330.     t[Z] = v[X]*mat[X][Z] + v[Y]*mat[Y][Z] + v[Z]*mat[Z][Z] + mat[W][Z];
  331.     t[W] = v[X]*mat[X][W] + v[Y]*mat[Y][W] + v[Z]*mat[Z][W] + mat[W][W];
  332.     vcopy4(t, vt);
  333. }
  334.  
  335. void mcopy(const Matrix m1, Matrix m2) {
  336.     register int i, j;
  337.  
  338.     for (i = X; i <= W; i++)
  339.     for (j = X; j <= W; j++)
  340.         m2[i][j] = m1[i][j];
  341. }
  342.  
  343. void mbuild_polar_rot(const float p1[2], Matrix m) {
  344.  
  345.     m[X][X] = cosf(p1[THETA]);
  346.     m[X][Y] = 0.0;
  347.     m[X][Z] = sinf(p1[THETA]);
  348.     m[X][W] = 0.0;
  349.  
  350.     m[Y][X] = sinf(p1[PHI]) * sinf(p1[THETA]);
  351.     m[Y][Y] = cosf(p1[PHI]);
  352.     m[Y][Z] = -sinf(p1[PHI]) * cosf(p1[THETA]);
  353.     m[Y][W] = 0.0;
  354.     
  355.     m[Z][X] = -cosf(p1[PHI]) * sinf(p1[THETA]);
  356.     m[Z][Y] = sinf(p1[THETA]);
  357.     m[Z][Z] = cosf(p1[PHI]) * cosf(p1[THETA]);
  358.     m[Z][W] = 0.0;
  359.     
  360.     m[W][X] = 0.0;
  361.     m[W][Y] = 0.0;
  362.     m[W][Z] = 0.0;
  363.     m[W][X] = 1.0;
  364. }
  365.  
  366. void mmake_rot_axis(const float a[3], const float alpha, Matrix m) {
  367.  
  368.     float q0, q1, q2, q3;
  369.     float s;
  370.     /* first build the quarternion */
  371.     s = sinf(alpha/2.0) / vlength(a);
  372.     q0 = a[X] * s;
  373.     q1 = a[Y] * s;
  374.     q2 = a[Z] * s;
  375.     q3 = cosf(alpha/2.0);
  376.     /* then build the matrix */
  377.     m[X][X] = 1.0 - 2.0 * (q1 * q1 + q2 * q2);
  378.     m[X][Y] = 2.0 * (q0 * q1 - q2 * q3);
  379.     m[X][Z] = 2.0 * (q2 * q0 + q1 * q3);
  380.     m[X][W] = 0.0;
  381.  
  382.     m[Y][X] = 2.0 * (q0 * q1 + q2 * q3);
  383.     m[Y][Y] = 1.0 - 2.0 * (q2 * q2 + q0 * q0);
  384.     m[Y][Z] = 2.0 * (q1 * q2 - q0 * q3);
  385.     m[Y][W] = 0.0;
  386.  
  387.     m[Z][X] = 2.0 * (q2 * q0 - q1 * q3);
  388.     m[Z][Y] = 2.0 * (q1 * q2 + q0 * q3);
  389.     m[Z][Z] = 1.0 - 2.0 * (q1 * q1 + q0 * q0);
  390.     m[Z][W] = 0.0;
  391.  
  392.     m[W][X] = 0.0;
  393.     m[W][Y] = 0.0;
  394.     m[W][Z] = 0.0;
  395.     m[W][W] = 1.0;
  396. }
  397.  
  398. void vrot_axis(const float v1[3], const float a[3], const float alpha,
  399.         float v2[3]) {
  400.     float q0, q1, q2, q3;
  401.     float s, tx, ty;
  402.     /* first build the quarternion */
  403.     s = sinf(alpha/2.0) / vlength(a);
  404.     q0 = a[X] * s;
  405.     q1 = a[Y] * s;
  406.     q2 = a[Z] * s;
  407.     q3 = cosf(alpha/2.0);
  408.     /* then calculate the new vector */
  409.     tx = (1.0 - 2.0 * (q1 * q1 + q2 * q2)) * v1[X] +
  410.         (2.0 * (q0 * q1 - q2 * q3)) * v1[Y] +
  411.         (2.0 * (q2 * q0 + q1 * q3)) * v1[Z];
  412.     ty = (2.0 * (q0 * q1 + q2 * q3)) * v1[X] +
  413.         (1.0 - 2.0 * (q2 * q2 + q0 * q0)) * v1[Y] +
  414.         (2.0 * (q1 * q2 - q0 * q3)) * v1[Z];
  415.     v2[Z] = (2.0 * (q2 * q0 - q1 * q3)) * v1[X] +
  416.         (2.0 * (q1 * q2 + q0 * q3)) * v1[Y] +
  417.         (1.0 - 2.0 * (q1 * q1 + q0 * q0)) * v1[Z];
  418.     v2[X] = tx;
  419.     v2[Y] = ty;
  420. }
  421.  
  422. /*
  423. *      rotate normalized vector a into normalized vector b.
  424. *      put the 4 by 4 rotation matrix in mat.
  425. *
  426. */
  427. void mmakerot(const float a[3], const float b[3], Matrix mat) {
  428.     float t1[4], t2[4], t3[4];
  429.     Matrix mat1, mat2;
  430.     int i;
  431.  
  432.     vcross(b, a, t1);
  433.     vnormal(t1);
  434.     vcross(b,t1,t2);
  435.     vcross(a,t1,t3);
  436.     vnormal(t2);
  437.     vnormal(t3);
  438.     mat1[W][W] = mat2[W][W] = 1.0;
  439.     for(i=X; i<=Z; i++) {
  440.     mat1[i][W] = mat1[W][i] = mat2[i][W] = mat2[W][i] = 0.0;
  441.     mat1[X][i] = b[i];
  442.     mat1[Y][i] = t2[i];
  443.     mat1[Z][i] = t1[i];
  444.     mat2[i][X] = a[i];
  445.     mat2[i][Y] = t3[i];
  446.     mat2[i][Z] = t1[i];
  447.     }
  448.     mmultmatrix(mat1,mat2,mat);
  449. }
  450.  
  451. void minvert(const Matrix mat, Matrix result) {
  452.     int i, j, k;
  453.     double temp;
  454.     double m[8][4];
  455.     /*   Declare identity matrix   */
  456.  
  457.     mcopy(idmatrix, result);
  458.     for (i = X;  i <= W;  i++) {
  459.     for (j = X;  j <= W;  j++) {
  460.         m[i][j] = mat[i][j];
  461.         m[i+4][j] = result[i][j];
  462.     }
  463.     }
  464.  
  465.     /*   Work across by columns   */
  466.  
  467.     for (i = X;  i <= W;  i++) {
  468.     for (j = i;  (m[i][j] == 0.0) && (j <= W);  j++) {
  469.         ;;            /* avoid warning empty for */
  470.     }
  471.     if (j == 4) {
  472.         fprintf (stderr, "error:  cannot do inverse matrix\n");
  473.         exit (2);
  474.     } 
  475.     else if (i != j) {
  476.         for (k = 0;  k < 8;  k++) {
  477.         temp = m[k][i];   
  478.         m[k][i] = m[k][j];   
  479.         m[k][j] = temp;
  480.         }
  481.     }
  482.  
  483.     /*   Divide original row   */
  484.  
  485.     for (j = 7;  j >= i;  j--)
  486.         m[j][i] /= m[i][i];
  487.  
  488.     /*   Subtract other rows   */
  489.  
  490.     for (j = X;  j <= W;  j++)
  491.         if (i != j)
  492.         for (k = 7;  k >= i;  k--)
  493.             m[k][j] -= m[k][i] * m[i][j];
  494.     }
  495.  
  496.     for (i = X;  i <= W;  i++)
  497.     for (j = X;  j <= W;  j++)
  498.         result[i][j] = m[i+4][j];
  499. }
  500.  
  501. /*
  502. * Get combined Model/View/Projection matrix, in any mmode
  503. */
  504. void mgetmatrix(Matrix m) {
  505.     long mm;
  506.  
  507.     mm = getmmode();
  508.  
  509.     if (mm == MSINGLE)
  510.     {
  511.     getmatrix(m);
  512.     }
  513.     else
  514.     {
  515.     Matrix mp, mv;
  516.  
  517.     mmode(MPROJECTION);
  518.     getmatrix(mp);
  519.     mmode(MVIEWING);
  520.     getmatrix(mv);
  521.  
  522.     pushmatrix();            /* Multiply them together */
  523.     loadmatrix(mp);
  524.     multmatrix(mv);
  525.     getmatrix(m);
  526.     popmatrix();
  527.  
  528.     mmode((short)mm);        /* Back into the mode we started in */
  529.     }
  530. }
  531.  
  532. /*
  533. * Gaussian Elimination with Scaled Column Pivoting
  534. *
  535. * copied out of the book by        Wade Olsen
  536. *                    Silicon Graphics
  537. *                    Feb. 12, 1990
  538. */
  539.  
  540. void linsolve(    
  541.     const float *eqs[],     /* System of eq's to solve */
  542.     int        n,         /* of size inmat[n][n+1] */
  543.     float    *x        /* Result float *or of size x[n] */
  544. )
  545. {
  546.     int        i, j, p;
  547.  
  548.     float **a;
  549.  
  550.     /* Allocate space to work in */
  551.     /* (avoid modifying the equations passed) */
  552.     a = (float **)malloc(sizeof(float *)*n);
  553.     for (i = 0; i < n; i++)
  554.     {
  555.     a[i] = (float *)malloc(sizeof(float)*(n+1));
  556.     bcopy(eqs[i], a[i], sizeof(float)*(n+1));
  557.     }
  558.  
  559.  
  560.     if (n == 1) {        /* The simple single variable case */
  561.     x[0] = a[0][1] / a[0][0];
  562.     return;
  563.     }        
  564.     /* Gausian elimination process */    
  565.     for (i = 0; i < n -1; i++) {
  566.     /* find non-zero pivot row */
  567.     p = i;
  568.     while (a[p][i] == 0.0) {
  569.         p++;
  570.         if (p == n) {
  571.         printf("linsolv:  No unique solution exists.\n");
  572.         exit(1);
  573.         }
  574.     }
  575.     /* row swap */
  576.     if (p != i) {
  577.         float *swap;
  578.  
  579.         swap = a[i];
  580.         a[i] = a[p];
  581.         a[p] = swap;
  582.     }
  583.     /* row subtractions */
  584.     for (j = i + 1; j < n; j++) {
  585.         float mult    = a[j][i] / a[i][i];
  586.  
  587.         int    k;
  588.         for (k = i + 1; k < n + 1; k++)
  589.         a[j][k] -= mult * a[i][k];
  590.     }
  591.     }
  592.  
  593.     if (a[n-1][n-1] == 0.0) {
  594.     printf("linsolv:  No unique solution exists.\n");
  595.     exit(1);
  596.     }
  597.  
  598.     /* backwards substitution */
  599.     x[n-1] = a[n-1][n] / a[n-1][n-1];
  600.     for (i = n -2; i > -1; i--) {
  601.     float sum = a[i][n];
  602.  
  603.     for (j = i + 1; j < n; j++)
  604.         sum -= a[i][j] * x[j];
  605.  
  606.     x[i] = sum / a[i][i];
  607.     }
  608.  
  609.     /* Free working space */
  610.     for (i = 0; i < n; i++)
  611.     {
  612.     free(a[i]);
  613.     }
  614.     free(a);
  615. }
  616.